
#include <iostream>
#include <stdlib.h>
#include <math.h>
#include <GL/glew.h>
#include <GL/glut.h>
#include <cuda_runtime.h>
#include <cutil.h>
#include <cuda_gl_interop.h>
#include <assert.h>
#include "Calculos.h"
#include "RadixSort/radixsort.h"
#include "Scan/Scan.h"
#include "SPHKernels.h"

using namespace std;


void radixSort() {
    int keybits = 32;
    RadixSort radixsort(totalParticulas, false);
    radixsort.sort(keysDevice, valuesDevice, totalParticulas, keybits);
}

void scan() {
    Scan scanner(totalParticulas);
    scanner.CriaOffset(countDevice, offSetDevice, totalCelulasGrid);
}

float calculaMassaEspecifica(int i, int j, float h, int indiceAtual) {
    float massaEsp = 0;
    int idCell = (i * resolucaoGridX) + j;

        int kInicial = offSet[idCell];
        int kFinal = offSet[idCell] + count[idCell];

        //cout << "\t\t i  = " << i << "\t j = " << j << "\tidCell = " << idCell
        //       << "\t kInicial = " << kInicial << "\t kFinal = " << kFinal << endl;

        for (int k = kInicial; k < kFinal; k++) {
            int indiceParticulas = values[k];
            if (indiceParticulas != indiceAtual) {
                //cout << "\t\t Vizinho (massa esp.) => " << indiceParticulas << endl;
                float distancia = sqrt((posicoesParticulasT[indiceAtual].x - posicoesParticulasT[indiceParticulas].x) * (posicoesParticulasT[indiceAtual].x - posicoesParticulasT[indiceParticulas].x) + (posicoesParticulasT[indiceAtual].y - posicoesParticulasT[indiceParticulas].y) * (posicoesParticulasT[indiceAtual].y - posicoesParticulasT[indiceParticulas].y));
                //cout << "\t\tDistancia para o ponto " << indiceParticulas << " foi: " << distancia << "\t h = " << h << endl;
                if (distancia < h) {
                    //cout << "\t\t (Distancia menor que h)" << endl;
                    massaEsp += massa * kernelPoly6(h, distancia);
                    //cout << "\t\tmassa especifica da particula " << indiceParticulas << " eh " << massa * (315 / (64 * 3.14 * pow(h, 9))) * pow((h * h - distancia * distancia), 3) << endl;
                }
            }
        }
    //cout << "\t\t Valor da massa esp. foi " << massaEsp << endl;
    return massaEsp;
}

float3 calculaFocaPressao(int i, int j, float h, int indiceAtual) {
    float3 forcaPressao;
    forcaPressao.x = 0;
    forcaPressao.y = 0;
    forcaPressao.z = 0.0f;
    int idCell = (i * resolucaoGridX) + j;
    float h2 = h*h;

    int kInicial = offSet[idCell];
    int kFinal = offSet[idCell] + count[idCell];

    //cout << "\t\t i  = " << i << "\t j = " << j << "\tidCell = " << idCell
    //     << "\t kInicial = " << kInicial << "\t kFinal = " << kFinal << endl;

    for (int k = kInicial; k < kFinal; k++) {
        int indiceParticulas = values[k];
        if (indiceParticulas != indiceAtual) {
            //cout << "\t\t Vizinho (pressao) => " << indiceParticulas << endl;
            float distancia2 = (posicoesParticulasT[indiceAtual].x - posicoesParticulasT[indiceParticulas].x) * (posicoesParticulasT[indiceAtual].x - posicoesParticulasT[indiceParticulas].x) + (posicoesParticulasT[indiceAtual].y - posicoesParticulasT[indiceParticulas].y) * (posicoesParticulasT[indiceAtual].y - posicoesParticulasT[indiceParticulas].y);
            //float distanciaX = (posicoesParticulasT[indiceAtual].x - posicoesParticulasT[indiceParticulas].x);
            //float distanciaY = (posicoesParticulasT[indiceAtual].y - posicoesParticulasT[indiceParticulas].y);
            //float distanciaZ = (posicoesParticulasT[indiceAtual].z - posicoesParticulasT[indiceParticulas].z);
            assert(distancia2 >= 0);
            if (distancia2 < h2) {
                float distancia = sqrt(distancia2);
                //cout << "\t\t Distancia(" << distancia << ") menor que h(" << h << ")" << endl;
                float3 gradientePressao;
                float mod;
                mod = sqrt(pow(posicoesParticulasT[indiceAtual].x - posicoesParticulasT[indiceParticulas].x, 2) +
                        pow(posicoesParticulasT[indiceAtual].y - posicoesParticulasT[indiceParticulas].y, 2) +
                        pow(posicoesParticulasT[indiceAtual].z - posicoesParticulasT[indiceParticulas].z, 2));
                gradientePressao.x = (posicoesParticulasT[indiceParticulas].x - posicoesParticulasT[indiceAtual].x)/mod;
                gradientePressao.y = (posicoesParticulasT[indiceParticulas].y - posicoesParticulasT[indiceAtual].y)/mod;
                gradientePressao.z = (posicoesParticulasT[indiceParticulas].z - posicoesParticulasT[indiceAtual].z)/mod;

                float kernelPress = kernelSpiky(h, distancia);
                forcaPressao.x += massa * ((particulasT2[indiceParticulas].pressao + particulasT[indiceAtual].pressao) / 2 * particulasT[indiceParticulas].massaEspecifica) * gradientePressao.x * kernelPress;
                //cout << "\t\t Pressao da particula " << indiceParticulas << " eh " << particulasT[indiceParticulas].pressao << "\t Massa esp. da particula " << indiceParticulas << " eh " << particulasT[indiceParticulas].massaEspecifica << endl;
                //cout << "\t\tforcaPressao.x da particula " << indiceParticulas << " eh " << forcaPressao.x << endl;
                forcaPressao.y += massa * ((particulasT2[indiceParticulas].pressao + particulasT[indiceAtual].pressao) / 2 * particulasT[indiceParticulas].massaEspecifica) * gradientePressao.y * kernelPress;
                //cout << "\t\tforcaPressao.y da particula " << indiceParticulas << " eh " << forcaPressao.y << endl;
                //forcaPressao.z += massa * ((particulasT[indiceParticulas].pressao + particulasT[indiceAtual].pressao) / 2 * particulasT[indiceParticulas].massaEspecifica) * (15 / 3.14 * pow(h, 6) * (pow((h - distanciaZ), 3)));
            }
        }
    }
    return forcaPressao;
}

float3 calculaFocaVisc(int i, int j, float h, int indiceAtual) {
    float3 forcaVisc;
    forcaVisc.x = 0;
    forcaVisc.y = 0;
    forcaVisc.z = 0;
    int idCell = (i * resolucaoGridX) + j;

    int kInicial = offSet[idCell];
    int kFinal = offSet[idCell] + count[idCell];

    //cout << "\t\t i  = " << i << "\t j = " << j << "\tidCell = " << idCell
    //     << "\t kInicial = " << kInicial << "\t kFinal = " << kFinal << endl;

    for (int k = kInicial; k < kFinal; k++) {
        int indiceParticulas = values[k];
        if (indiceParticulas != indiceAtual) {
            //cout << "\t\t Vizinho (visc.) => " << indiceParticulas << endl;
            float distancia = sqrt((posicoesParticulasT[indiceAtual].x - posicoesParticulasT[indiceParticulas].x) * (posicoesParticulasT[indiceAtual].x - posicoesParticulasT[indiceParticulas].x) + (posicoesParticulasT[indiceAtual].y - posicoesParticulasT[indiceParticulas].y) * (posicoesParticulasT[indiceAtual].y - posicoesParticulasT[indiceParticulas].y));
            if (distancia > 0 && distancia < h) {
                //cout << "\t\t (Distancia menor que h)" << endl;
                forcaVisc.x += massa * ((particulasT[indiceParticulas].velocidade.x + particulasT[indiceAtual].velocidade.x) / 2 * particulasT[indiceParticulas].massaEspecifica) * kernelViscosity(h, distancia);
                forcaVisc.y += massa * ((particulasT[indiceParticulas].velocidade.y + particulasT[indiceAtual].velocidade.y) / 2 * particulasT[indiceParticulas].massaEspecifica) * kernelViscosity(h, distancia);
                //forcaVisc.z += massa * ((particulasT[indiceParticulas].velocidade.z + particulasT[indiceAtual].velocidade.z) / 2 * particulasT[indiceParticulas].massaEspecifica) * (45/3.14*pow(h,6)* h - distanciaZ);
            }
        }
    }
    //cout << "\t\t Valor da viscosidade.x foi " << forcaVisc.x << "\t viscosidade.y foi " << forcaVisc.y << "\t viscosidade.z foi " << forcaVisc.z << endl;
    return forcaVisc;
}

void corrigePosicao(int posicaoParticulaAtual){

    if(posicoesParticulasT2[posicaoParticulaAtual].x < XminGrid){
        posicoesParticulasT2[posicaoParticulaAtual].x = XminGrid + ABS(posicoesParticulasT2[posicaoParticulaAtual].x - XminGrid);
    }
    if(posicoesParticulasT2[posicaoParticulaAtual].x > XmaxGrid){
        posicoesParticulasT2[posicaoParticulaAtual].x = XmaxGrid - ABS(posicoesParticulasT2[posicaoParticulaAtual].x - XmaxGrid);
    }
    if(posicoesParticulasT2[posicaoParticulaAtual].y < YminGrid){
        posicoesParticulasT2[posicaoParticulaAtual].y = YminGrid + ABS(posicoesParticulasT2[posicaoParticulaAtual].y - YminGrid);
    }
    if(posicoesParticulasT2[posicaoParticulaAtual].y > YmaxGrid){
        posicoesParticulasT2[posicaoParticulaAtual].y = YmaxGrid + ABS(posicoesParticulasT2[posicaoParticulaAtual].y - YmaxGrid);
    }

    /*if(posicoesParticulasT2[posicaoParticulaAtual].x < (XminGrid + recuoXGrid)){
        float novaPosicao = ABS((XminGrid + recuoXGrid) - posicoesParticulasT2[posicaoParticulaAtual].x);
        if(novaPosicao > ((XmaxGrid - recuoXGrid) - (XminGrid + recuoXGrid))){
            novaPosicao = (XmaxGrid - recuoXGrid) - (XminGrid + recuoXGrid);
        }
        //posicoesParticulasT2[posicaoParticulaAtual].x = (XminGrid + recuoXGrid) + novaPosicao;
        posicoesParticulasT2[posicaoParticulaAtual].x = (XmaxGrid - recuoXGrid) - novaPosicao;
        //posicoesParticulasT2[posicaoParticulaAtual].x = (XmaxGrid - recuoXGrid - recuoXGrid);

    }
    if(posicoesParticulasT2[posicaoParticulaAtual].x > (XmaxGrid - recuoXGrid)){
        float novaPosicao = ABS(posicoesParticulasT2[posicaoParticulaAtual].x - (XmaxGrid - recuoXGrid));
        if(novaPosicao > ((XmaxGrid - recuoXGrid) - (XminGrid + recuoXGrid))){
            novaPosicao = (XmaxGrid - recuoXGrid) - (XminGrid + recuoXGrid);
        }
        //posicoesParticulasT2[posicaoParticulaAtual].x = (XmaxGrid - recuoXGrid) - novaPosicao;
        posicoesParticulasT2[posicaoParticulaAtual].x = (XminGrid + recuoXGrid) + novaPosicao;
        //posicoesParticulasT2[posicaoParticulaAtual].x = (XminGrid + recuoXGrid + recuoXGrid);
    }
    if(posicoesParticulasT2[posicaoParticulaAtual].y < (YminGrid + recuoYGrid)){
        float novaPosicao = ABS((YminGrid + recuoYGrid) - posicoesParticulasT2[posicaoParticulaAtual].y);
        if(novaPosicao > ((YmaxGrid - recuoYGrid) - (YminGrid + recuoYGrid))){
            novaPosicao = (YmaxGrid - recuoYGrid) - (YminGrid + recuoYGrid);
        }
        //posicoesParticulasT2[posicaoParticulaAtual].y = (YminGrid + recuoYGrid) + novaPosicao;
        posicoesParticulasT2[posicaoParticulaAtual].y = (YmaxGrid - recuoYGrid) - novaPosicao;
        //posicoesParticulasT2[posicaoParticulaAtual].y = (YmaxGrid - recuoYGrid - recuoYGrid);
    }
    if(posicoesParticulasT2[posicaoParticulaAtual].y > (YmaxGrid - recuoYGrid)){
        float novaPosicao = ABS(posicoesParticulasT2[posicaoParticulaAtual].y - (YmaxGrid - recuoYGrid));
        if(novaPosicao > ((YmaxGrid - recuoYGrid) - (YminGrid + recuoYGrid))){
            novaPosicao = (YmaxGrid - recuoYGrid) - (YminGrid + recuoYGrid);
        }
        //posicoesParticulasT2[posicaoParticulaAtual].y = (YmaxGrid - recuoYGrid) - novaPosicao;
        posicoesParticulasT2[posicaoParticulaAtual].y = (YminGrid + recuoYGrid) + novaPosicao;
        //posicoesParticulasT2[posicaoParticulaAtual].y = (YminGrid + recuoYGrid + recuoYGrid);
    }*/

}

void defineVetoresTeT2(){
    if(VBO_Ativo){
        posicoesParticulasT = posicoesParticulas;
        posicoesParticulasT2 = posicoesParticulas2;
        particulasT = particulas;
        particulasT2 = particulas2;
    } else {
        posicoesParticulasT = posicoesParticulas2;
        posicoesParticulasT2 = posicoesParticulas;
        particulasT = particulas2;
        particulasT2 = particulas;
    }
}

void calculaParte1() {
    if (VBO_Ativo) {
        // Map buffer object for writing from CUDA
        CUDA_SAFE_CALL(cudaGLMapBufferObject((void**) & VBO, positionsVBO));
    } else {
        // Map buffer object for writing from CUDA
        CUDA_SAFE_CALL(cudaGLMapBufferObject((void**) & VBO2, positionsVBO2));
    }

    launch_ZeraCountEOffSet(blocksPerGrid, threadsPerBlock, totalCelulasGrid, countDevice, offSetDevice);

    if (VBO_Ativo) {
        launch_CalculaIDGridCell(blocksPerGrid, threadsPerBlock, VBO, valuesDevice,
                keysDevice, totalParticulas, resolucaoGridX, deltaXGrid, deltaYGrid, countDevice,
                XminGrid, YminGrid);
    } else {
        launch_CalculaIDGridCell(blocksPerGrid, threadsPerBlock, VBO2, valuesDevice,
                keysDevice, totalParticulas, resolucaoGridX, deltaXGrid, deltaYGrid, countDevice,
                XminGrid, YminGrid);
    }

    scan();

    // Copy result from device memory to host memory
    CUDA_SAFE_CALL(cudaMemcpy(keys, keysDevice, totalParticulas * sizeof (unsigned int), cudaMemcpyDeviceToHost));
    CUDA_SAFE_CALL(cudaMemcpy(count, countDevice, totalCelulasGrid * sizeof (unsigned int), cudaMemcpyDeviceToHost));
    CUDA_SAFE_CALL(cudaMemcpy(offSet, offSetDevice, totalCelulasGrid * sizeof (unsigned int), cudaMemcpyDeviceToHost));

    radixSort();
    CUDA_SAFE_CALL(cudaMemcpy(keys, keysDevice, totalParticulas * sizeof (unsigned int), cudaMemcpyDeviceToHost));
    CUDA_SAFE_CALL(cudaMemcpy(values, valuesDevice, totalParticulas * sizeof (unsigned int), cudaMemcpyDeviceToHost));
    CUDA_SAFE_CALL(cudaMemcpy(particulas, particulasDevice, totalParticulas * sizeof (Particula), cudaMemcpyDeviceToHost));
    CUDA_SAFE_CALL(cudaMemcpy(particulas2, particulasDevice2, totalParticulas * sizeof (Particula), cudaMemcpyDeviceToHost));

    if (VBO_Ativo) {
        // Unmap buffer object
        CUDA_SAFE_CALL(cudaGLUnmapBufferObject(positionsVBO));
    } else {
        // Unmap buffer object
        CUDA_SAFE_CALL(cudaGLUnmapBufferObject(positionsVBO2));
    }
}

void calculaOnCPU(void) {
    defineVetoresTeT2();
    calculaParte1();
    deleteVBO();
    //glBindBufferARB(GL_ARRAY_BUFFER_ARB, positionsVBO);
    //float3* ptr = (float3*) glMapBufferARB(GL_ARRAY_BUFFER_ARB, GL_READ_WRITE_ARB);

    for (int indiceValue = 0; indiceValue < totalParticulas; indiceValue++) {
        int jAtual = 0;
        int iAtual = 0;
        int jMin = 0;
        int jMax = resolucaoGridX;
        int iMin = 0;
        int iMax = resolucaoGridY;
        float h = deltaXGrid * 1.1; //sqrt(pow(deltaXGrid, 2) + pow(deltaYGrid, 2));
        //float h = deltaXPart * 1.502;
        //if (deltaXPart > deltaYPart) {
        //    h = deltaYPart * 1.502;
        //}
        cout << "h = " << h << "\t DeltaXpart = " << deltaXPart << "\t DeltaYpart = " << deltaYPart << endl;
        float3 forcaPressao = {0.0f, 0.0f, 0.0f};
        float3 forcaAceleracao = {0.0f, 0.0f, 0.0f};
        float3 forcaViscosidade = {0.0f, 0.0f, 0.0f};
        float3 gravidade = {0.0f, -9.8f, 0.0f};
        /*forcaPressaoVetor[indiceValue].x = forcaPressaoVetor[indiceValue].y = forcaPressaoVetor[indiceValue].z = 0;
        forcaAceleracaoVetor[indiceValue].x = forcaAceleracaoVetor[indiceValue].y = forcaAceleracaoVetor[indiceValue].z = 0;
        forcaViscosidadeVetor[indiceValue].x = forcaViscosidadeVetor[indiceValue].y = forcaViscosidadeVetor[indiceValue].z = 0;*/

        cout << "\nPara o indice indice value igual a = " << indiceValue << " e values[indiceValue] = " << values[indiceValue] << endl;

        if (indiceValue < totalParticulas) {
            int posicaoParticulaAtual = values[indiceValue];
            jAtual = (posicoesParticulasT[posicaoParticulaAtual].x / deltaXGrid);
            iAtual = (posicoesParticulasT[posicaoParticulaAtual].y / deltaYGrid);
            //int idCellAtual = (iAtual * numeroColunasGridCell) + jAtual;
            cout << "\t jAtual = " << jAtual << "    iAtual = " << iAtual;

            //VERIFICAO DAS EXTREMIDADES.
            if (jAtual == 0) {
                jMin = 0;
            } else {
                jMin = jAtual - 1;
            }

            if (jAtual == (resolucaoGridX - 1)) {
                jMax = jAtual;
            } else {
                jMax = jAtual + 1;
            }

            if (iAtual == 0) {
                iMin = 0;
            } else {
                iMin = iAtual - 1;
            }

            if (iAtual == (resolucaoPartY - 1)) {
                iMax = iAtual;
            } else {
                iMax = iAtual + 1;
            }

            cout << "\t (jMin = " << jMin << "    jMax = " << jMax << "    iMin = " << iMin << "    iMax = " << iMax << ")" << endl;
            particulasT2[posicaoParticulaAtual].massaEspecifica = 0;

            //CALCULA AS MASSAS ESPECIFICAS.
            for (int i = iMin; i <= iMax; i++) {
                for (int j = jMin; j <= jMax; j++) {
                    particulasT2[posicaoParticulaAtual].massaEspecifica += calculaMassaEspecifica(i, j, h, posicaoParticulaAtual);
                    //cout << "\t\t Massa esp. temp da particula " << posicaoParticulaAtual << " eh " << particulasT1[posicaoParticulaAtual].massaEspecifica << endl;
                }
            }

            cout << "\t pos.x: " << posicoesParticulasT[posicaoParticulaAtual].x
                    << "\t pos.y: " << posicoesParticulasT[posicaoParticulaAtual].y
                    << "\t pos.z: " << posicoesParticulasT[posicaoParticulaAtual].z << endl;
            cout << "\t Massa especifica final da particula " << posicaoParticulaAtual << " eh " << particulasT2[posicaoParticulaAtual].massaEspecifica << endl ;

            //CALCALA A PRESSAO NA PARTICULA.
            particulasT2[posicaoParticulaAtual].pressao = constGas * (particulasT2[posicaoParticulaAtual].massaEspecifica - massaEspecificaInicial);

            //CALCULA A FORÇA DE PRESSAO.
            for (int i = iMin; i <= iMax; i++) {
                for (int j = jMin; j <= jMax; j++) {
                    float3 forcaPressaoTemp = calculaFocaPressao(i, j, h, posicaoParticulaAtual);
                    forcaPressao.x += forcaPressaoTemp.x;
                    forcaPressao.y += forcaPressaoTemp.y;
                    forcaPressao.z += forcaPressaoTemp.z;
                }
            }
            forcaPressao.x *= -1;
            forcaPressao.y *= -1;
            forcaPressao.z *= -1;
            forcaPressaoVetor[posicaoParticulaAtual].x = forcaPressao.x;
            forcaPressaoVetor[posicaoParticulaAtual].y = forcaPressao.y;
            forcaPressaoVetor[posicaoParticulaAtual].z = forcaPressao.z;
            cout << "\t ForcaPressao.x: " << forcaPressao.x << "\t ForcaPressao.y: " << forcaPressao.y << "\t ForcaPressao.z: " << forcaPressao.z << endl;

            //CALCULA A FORÇA VISCOSIDADE.
            for (int i = iMin; i <= iMax; i++) {
                for (int j = jMin; j <= jMax; j++) {
                    float3 forcaViscTemp = calculaFocaVisc(i, j, h, posicaoParticulaAtual);
                    forcaViscosidade.x += forcaViscTemp.x;
                    forcaViscosidade.y += forcaViscTemp.y;
                    forcaViscosidade.z += forcaViscTemp.z;
                }
            }
            forcaViscosidade.x *= coeficienteViscosidade;
            forcaViscosidade.y *= coeficienteViscosidade;
            forcaViscosidade.z *= coeficienteViscosidade;
            forcaViscosidadeVetor[posicaoParticulaAtual].x = forcaViscosidade.x;
            forcaViscosidadeVetor[posicaoParticulaAtual].y = forcaViscosidade.y;
            forcaViscosidadeVetor[posicaoParticulaAtual].z = forcaViscosidade.z;
            cout << "\t ForcaViscosidade.x: " << forcaViscosidade.x << "\t ForcaViscosidade.y: " << forcaViscosidade.y << "\t ForcaViscosidade.z: " << forcaViscosidade.z << endl;

            forcaAceleracao.x = forcaPressao.x + forcaViscosidade.x; //+ particulasT1[posicaoParticulaAtual].massaEspecifica * gravidade.x;
            forcaAceleracao.y = forcaPressao.y + forcaViscosidade.y; //+ particulasT1[posicaoParticulaAtual].massaEspecifica * gravidade.y;
            forcaAceleracao.z = forcaPressao.z + forcaViscosidade.z; //+ particulasT1[posicaoParticulaAtual].massaEspecifica * gravidade.z;
            forcaAceleracaoVetor[posicaoParticulaAtual].x = forcaAceleracao.x;
            forcaAceleracaoVetor[posicaoParticulaAtual].y = forcaAceleracao.y;
            forcaAceleracaoVetor[posicaoParticulaAtual].z = forcaAceleracao.z;
            cout << "\t ForcaAceleracao.x: " << forcaAceleracao.x << "\t ForcaAceleracao.y: " << forcaAceleracao.y << "\t ForcaAceleracao.z: " << forcaAceleracao.z << endl;

            particulasT2[posicaoParticulaAtual].velocidade.x = particulasT[posicaoParticulaAtual].velocidade.x + passoIntegracao / 2 * forcaAceleracao.x;
            particulasT2[posicaoParticulaAtual].velocidade.y = particulasT[posicaoParticulaAtual].velocidade.y + passoIntegracao / 2 * forcaAceleracao.y;
            particulasT2[posicaoParticulaAtual].velocidade.z = 0;//particulasT[posicaoParticulaAtual].velocidade.z + passoIntegracao / 2 * forcaAceleracao.z;
            cout << "\t Velocidade.x: " << particulasT2[posicaoParticulaAtual].velocidade.x
                    << "\t Velocidade.y: " << particulasT2[posicaoParticulaAtual].velocidade.y
                    << "\t Velocidade.z: " << particulasT2[posicaoParticulaAtual].velocidade.z << endl;

            posicoesParticulasT2[posicaoParticulaAtual].x = posicoesParticulasT[posicaoParticulaAtual].x + passoIntegracao / 2 * particulasT2[posicaoParticulaAtual].velocidade.x;
            cout << passoIntegracao / 2 * particulasT2[posicaoParticulaAtual].velocidade.x;
            posicoesParticulasT2[posicaoParticulaAtual].y = posicoesParticulasT[posicaoParticulaAtual].y + passoIntegracao / 2 * particulasT2[posicaoParticulaAtual].velocidade.y;
            cout << "\t" << passoIntegracao / 2 * particulasT2[posicaoParticulaAtual].velocidade.x << endl;
            posicoesParticulasT2[posicaoParticulaAtual].z = 0;//posicoesParticulasT[posicaoParticulaAtual].z + passoIntegracao / 2 * particulasT1[posicaoParticulaAtual].velocidade.z;

            corrigePosicao(posicaoParticulaAtual);

            cout << "\t pos.x: " << posicoesParticulasT2[posicaoParticulaAtual].x
                    << "\t pos.y: " << posicoesParticulasT2[posicaoParticulaAtual].y
                    << "\t pos.z: " << posicoesParticulasT2[posicaoParticulaAtual].z << endl;

            /*int jTemp = (posicoesParticulasT2[posicaoParticulaAtual].x) / deltaXGrid;
            int iTemp = (posicoesParticulasT2[posicaoParticulaAtual].y) / deltaYGrid;
            int idCell = (iTemp * numeroColunasGridCell) + jTemp;
            cout << indiceValue << "-> Next idCell: " << idCell << endl;*/

        }
    }

    if (VBO_Ativo) {
        VBO_Ativo = false;
    } else {
        VBO_Ativo = true;
    }

    buildVBO();

    glutPostRedisplay();
}

void CalculaKernelsOnGPU() {
    /*if (VBO_Ativo) {
        // Map buffer object for writing from CUDA
        CUDA_SAFE_CALL(cudaGLMapBufferObject((void**) & VBO, positionsVBO));

        launch_ZeraCountEOffSet(blocksPerGrid, threadsPerBlock, numeroCelulasGrid, countDevice, offSetDevice);

        launch_CalculaIDGridCell(blocksPerGrid, threadsPerBlock, VBO, valuesDevice,
                keysDevice, numeroParticulas, numeroColunasGridCell, deltaXGrid, deltaYGrid, countDevice,
                XminGrid, YminGrid);

        scan();

        // Copy result from device memory to host memory
        CUDA_SAFE_CALL(cudaMemcpy(keys, keysDevice, numeroParticulas * sizeof (unsigned int), cudaMemcpyDeviceToHost));
        //cudaMemcpy(values, valuesDevice, numeroParticulas*sizeof(unsigned int), cudaMemcpyDeviceToHost);
        CUDA_SAFE_CALL(cudaMemcpy(count, countDevice, numeroCelulasGrid * sizeof (unsigned int), cudaMemcpyDeviceToHost));
        CUDA_SAFE_CALL(cudaMemcpy(offSet, offSetDevice, numeroCelulasGrid * sizeof (unsigned int), cudaMemcpyDeviceToHost));

#ifdef DEBUG
        printCells();

        printCount();

        cout << "\n Fazendo o radixSort" << endl;
#endif

        radixSort();
        CUDA_SAFE_CALL(cudaMemcpy(keys, keysDevice, numeroParticulas * sizeof (unsigned int), cudaMemcpyDeviceToHost));
        CUDA_SAFE_CALL(cudaMemcpy(values, valuesDevice, numeroParticulas * sizeof (unsigned int), cudaMemcpyDeviceToHost));

#ifdef DEBUG
        for (int i = 0; i < numeroParticulas; i++) {
            if (i < numeroCelulasGrid) {
                cout << "id: " << values[i] << "\t idCell: " << keys[i] << "\t offSet: " << offSet[i]
                        << "\t count[" << i << "] : " << count[i] << endl;
            } else {
                cout << "id: " << values[i] << "\t idCell: " << keys[i] << endl;
            }
        }
//        getchar();
#endif

        launch_FindNeighbors(blocksPerGrid, threadsPerBlock, VBO, particulasDevice, particulasDevice2, valuesDevice, keysDevice,
                countDevice, offSetDevice, numeroColunasGridCell, numeroLinhasGridCell, numeroParticulas,
                deltaXGrid, deltaYGrid, XminGrid, YminGrid, massaEspecificaInicial, massa);

        //CUDA_SAFE_CALL(cudaMemcpy(particulas, particulasDevice, numeroParticulas * sizeof (Particula), cudaMemcpyDeviceToHost));
        CUDA_SAFE_CALL(cudaMemcpy(particulas2, particulasDevice2, numeroParticulas * sizeof (Particula), cudaMemcpyDeviceToHost));

#ifdef DEBUG
        cout << "Imprimindo particulas t (particulas)." << endl;
        printDadosParticulas();
        cout << "\n Imprimindo particulas t + 1 (particulas2)." << endl;
        printDadosParticulas2();
#endif
        // Unmap buffer object
        CUDA_SAFE_CALL(cudaGLUnmapBufferObject(positionsVBO));

        VBO_Ativo = false;

    } else {

        // Map buffer object for writing from CUDA
        CUDA_SAFE_CALL(cudaGLMapBufferObject((void**) & VBO, positionsVBO));

        launch_ZeraCountEOffSet(blocksPerGrid, threadsPerBlock, numeroCelulasGrid, countDevice, offSetDevice);

        launch_CalculaIDGridCell(blocksPerGrid, threadsPerBlock, VBO, valuesDevice,
                keysDevice, numeroParticulas, numeroColunasGridCell, deltaXGrid, deltaYGrid, countDevice,
                XminGrid, YminGrid);

        scan();

        // Copy result from device memory to host memory
        CUDA_SAFE_CALL(cudaMemcpy(keys, keysDevice, numeroParticulas * sizeof (unsigned int), cudaMemcpyDeviceToHost));
        //cudaMemcpy(values, valuesDevice, numeroParticulas*sizeof(unsigned int), cudaMemcpyDeviceToHost);
        CUDA_SAFE_CALL(cudaMemcpy(count, countDevice, numeroCelulasGrid * sizeof (unsigned int), cudaMemcpyDeviceToHost));
        CUDA_SAFE_CALL(cudaMemcpy(offSet, offSetDevice, numeroCelulasGrid * sizeof (unsigned int), cudaMemcpyDeviceToHost));

#ifdef DEBUG
        printCells();

        printCount();

        cout << "\n Fazendo o radixSort" << endl;
#endif

        radixSort();
        CUDA_SAFE_CALL(cudaMemcpy(keys, keysDevice, numeroParticulas * sizeof (unsigned int), cudaMemcpyDeviceToHost));
        CUDA_SAFE_CALL(cudaMemcpy(values, valuesDevice, numeroParticulas * sizeof (unsigned int), cudaMemcpyDeviceToHost));

#ifdef DEBUG
        for (int i = 0; i < numeroParticulas; i++) {
            if (i < numeroCelulasGrid) {
                cout << "id: " << values[i] << "\t idCell: " << keys[i] << "\t offSet: " << offSet[i]
                        << "\t count[" << i << "] : " << count[i] << endl;
            } else {
                cout << "id: " << values[i] << "\t idCell: " << keys[i] << endl;
            }
        }
//        getchar();
#endif

        launch_FindNeighbors(blocksPerGrid, threadsPerBlock, VBO, particulasDevice2, particulasDevice, valuesDevice, keysDevice,
                countDevice, offSetDevice, numeroColunasGridCell, numeroLinhasGridCell, numeroParticulas,
                deltaXGrid, deltaYGrid, XminGrid, YminGrid, massaEspecificaInicial, massa);

        //CUDA_SAFE_CALL(cudaMemcpy(particulas, particulasDevice, numeroParticulas * sizeof (Particula), cudaMemcpyDeviceToHost));
        CUDA_SAFE_CALL(cudaMemcpy(particulas, particulasDevice, numeroParticulas * sizeof (Particula), cudaMemcpyDeviceToHost));

#ifdef DEBUG
        cout << "Imprimindo particulas t + 1 (particulas)." << endl;
        printDadosParticulas();
        cout << "\n Imprimindo particulas t (particulas2)." << endl;
        printDadosParticulas2();
#endif
        // Unmap buffer object
        CUDA_SAFE_CALL(cudaGLUnmapBufferObject(positionsVBO));

        VBO_Ativo = true;
    }*/
}
